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1.  Project  Objectives 


•  Development  of  computationally  efficient  multiscale  modeling  techniques  for 
characterizing  the  damage  state  of  a  material  (including  nucleation  and  growth): 

o  Develop  physics  based  stochastic  framework  to  investigate  fatigue  damage  initiation 
and  growth 

o  Evaluate  potential  sites  for  damage  nucleation  in  the  microstructure 
o  Characterized  damage  evolution  on  complex  structures 

o  Formulate  methodologies  to  characterize  damage  evolution  at  different  length  scales 
and  stress  states 

o  Quantify  variability  of  microstructural  parameters  that  control  damage 
o  Incorporate  multiscale  modeling  with  virtual  sensing  techniques 
o  Improve  simulation  efficiency  and  accuracy  for  online  monitoring 
o  Incorporate  thermal  creep  effects  to  solve  variable  temperature  issues 

•  Investigation  of  damage  detection  and  classification  techniques  for  sensor  integration  and 
instrumentation: 

o  Identify  physics-based  nondestructive  methods  to  quantify  damage  on  the  microscopic 
and  mesoscopic  length  scales 

o  Perform  systems-level  analysis  using  both  sensor  data  and  physics-based  models 
o  Develop  sensor  management  schemes 

o  Design  optimal  detection  and  classification  algorithms  to  analyze,  classify  and  detect 
structural  damage 

o  Adapt  detection  and  classification  algorithms  to  ensure  robustness  to  material  or 
environment  variability 

o  Correlate  point-source/point-receiver  measurements  to  damage-induced  isotropy 
changes  in  selected  materials  systems 

o  Interpret  ultrasonic  results  within  broader  framework  of  orientation  distribution 
functions  and  crack-induced  texture 

o  Quantify  effects  of  temperature  on  damage  behavior  and  damage  monitoring  strategies 

•  Development  of  an  integrated  structural  health  monitoring  (SHM)  and  stochastic  prognosis 
framework  to  predict  failure  and  remaining  useful  life  of  metallic  components: 

o  Develop  stochastic  models  to  account  for  uncertainties  in  signals 
o  Construct  probabilistic  models  sensitive  to  fine  scale  information 
o  Use  modeling  information  for  optimal  sensing  networks 
o  Implement  models  for  estimating  residual  life 
o  Construct  quantitative  measures  of  damage 

o  Develop  stochastic  models  to  account  for  uncertainties  in  excitation  and  measured 
signals 

•  Testing,  validation  and  application: 

o  Characterize  damage  evolution  and  link  damage  parameters  to  sensor  output 
o  Validate  and  calibrate  modeling  techniques 

o  Characterize  damage  evolution  on  complex  sensor  integrated  structures 
o  Apply  methodology  to  structural  hot  spots 
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2.  Approach  and  Accomplishments 

Task  1:  Physics  based  Multiscale  Modeling 

Damage  assessment  and  residual  useful  life  estimation  (RULE)  are  essential  for 
aerospace,  civil  and  naval  structures.  Structural  Health  Monitoring  (SHM)  attempts  to  automate 
the  process  of  damage  detection  and  identification.  Multiscale  modeling  is  a  key  element  in 
SHM.  It  not  only  provides  important  information  on  the  physics  of  failure,  such  as  damage 
initiation  and  growth,  the  output  can  be  used  as  “virtual  sensing”  data  for  detection  and 
prognosis.  A  multiscale  model  has  been  developed  by  bridging  the  relevant  length  scales, 
starting  from  micro  to  meso  and  macro  (or  structural)  scale.  Micro  structural  representations 
obtained  from  material  characterization  are  used  to  define  the  length  scales  and  to  capture  the 
size  and  orientation  of  the  grains  at  the  micro  level.  A  microvoid  model  accounting  for  size  and 
crystal  orientation  effects  is  developed.  Parametric  studies  are  conducted  to  estimate  material 
parameters  used  in  this  constitutive  model.  Numerical  and  experimental  simulations  are 
performed  to  investigate  the  effects  of  Representative  Volume  Element  (RVE)  size,  defect  area 
fraction  and  distribution.  A  Statistical  Volume  Element  (SVE)  with  high  computational 
efficiency  is  also  developed  at  the  mesoscale  to  represent  the  microstructure  of  the  material.  A 
multiscale  damage  criterion  accounting  for  crystal  orientation  effect  is  developed  and  applied  for 
fatigue  crack  initial  stage  prediction.  The  micro  scale  damage  information  is  passed  from  local  to 
grain  level,  bridging  microscale  and  mesoscale,  in  the  form  of  a  damage  vector  via  averaging 
techniques.  The  damage  evolution  rule  for  a  meso  RVE,  which  contains  several  grains,  is 
calculated  by  modifying  the  Kreisselmeier-Steinhauser  (KS)  function,  which  is  used  in  multi- 
objective  optimization  applications,  to  obtain  lower  and  upper  bounds  of  damage  envelopes.  A 
critical  damage  value  is  derived  to  complete  the  damage  criterion  for  fatigue  life  prediction  in 
Aluminum  2024  test  articles. 

Material  characterization  and  testing 

Various  mechanical  testing  experiments  were  conducted  for  the  purpose  of  model 
calibration  and  validation.  Experimental  techniques  for  mechanical  testing  followed  ASTM 
standards  as  close  as  possible  depending  on  the  test  articles.  In  this  case,  tensile  testing  to 
determine  basic  mechanical  properties  under  monotonic  loading  followed  ASTM  E8.  The  results 
of  these  experiments,  for  dogbone  samples  machined  with  load  axes  either  parallel  or 
perpendicular  to  the  rolling  direction  of  the  A1  plate  used,  were  combined  with  hysteresis  loops 
obtained  under  load  control  for  similar  samples  tested  following  ASTM  E466  and  E606.  Fatigue 
crack  propagation  experiments  were  conducted  on  CT  specimens  following  ASTM  E647  to 
measure  crack  growth  kinetics  for  A1  and  Ti  alloy  samples  with  cracks  propagating  either 
parallel  or  perpendicular  to  the  rolling  direction  of  the  plate.  Testing  on  specimens  with  non¬ 
standard  geometries,  i.e.,  large  and  small  lug-joints,  as  well  as  large  and  small  cruciforms,  or 
under  nonstandard  testing  conditions,  e.g.,  constant  AK  experiments  to  propagate  fatigue  cracks 
in  CT  specimens  in  order  to  perform  in-situ  Digital  Image  Correlation  (DIC)  experiments,  were 
performed  following  best  practices  based  on  the  standards  mentioned  above  as  well  as 
experience.  Examples  of  the  geometry  of  some  of  the  samples  used  are  shown  in  Figure  1. 

Multiscale  characterization  of  damage  evolution  was  undertaken  using  a  dual  approach. 
First,  long  crack  propagation  was  characterized  to  study  correlations  with  microstructure  and 
applied  load.  Then,  crack  nucleation  and  short  crack  growth  were  studied  up  to  crack  lengths 
where  the  transition  to  long  cracks  occurs.  There  was  particular  interest  in  finding  correlations 
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Figure  1.  Typical  geometries  used  in  the  mechanical  testing,  (a)  CT  specimen  (W=25.4  mm), 
(b)  Notched  dogbone  (dimensions  in  inches),  (c)  Miniature  cruciform  sample  (sample  arms 
were  three  inches  long),  (d)  Large  lug-joint  (dimensions  in  inches).  The  small  lug-joint  had  all 
dimensions  scaled  down  by  a  factor  of  three. 


between  local  fields  ahead  or  around  defects  and  the  evolution  of  fatigue  damage.  In  the  case  of 
long  cracks,  data  analysis  techniques  developed  during  this  project  were  used  on  results  obtained 
elsewhere  for  strain  fields  around  fatigue  crack  tips  in  pure  polycrystalline  nickel.  The  results 
indicated  that  the  area  integral  of  the  opening  strain  ahead  of  the  crack  tip,  the  integrated  strain 
Sint,  could  be  used  as  non-local  strain  parameter  to  correlate  the  fields  ahead  of  crack  tip  with  the 
macroscopic  crack  growth  kinetics.  There  was  a  power-law  relationship  between  Sint  and  AK, 
where  the  exponent  was  equal  to  the  Paris  exponent  of  the  crack  growth  kinetics,  which  meant 
that  the  integrated  strain  was  proportional  to  the  crack  growth  rate.  Analysis  revealed  that  the 
Paris  exponent  m  could  be  predicted  directly  from  the  cyclic  hardening  exponent  n’  via  the 
following  relationship: 

m  =  2(i+Ii7)  (1) 

where  0  <  n’  <  1 .  This  gives  3  <  m  <  4,  which  accounts  for  about  50%  of  the  observed  values  of 
between  1.5  and  6  for  a  wide  variety  of  metallic  materials.  This  result  was  obtained  based  on  the 
fact  that  strain  ahead  of  the  fatigue  cracks  localized  along  well  defined  slip  bands  in  the 
equiaxed  grains  of  the  polycrystalline  Ni  samples  used,  where  the  integrated  strain  was 
proportional  to  the  accumulated  displacement  in  these  bands.  These  results  are  helpful  because 
they  offer  an  avenue  to  predict  crack  growth  kinetics  based  on  strain  fields  around  crack  tips. 
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which  can  be  calculated  using  mature  methodologies  like  finite  elements.  Prediction  of  crack 
growth  kinetics  is  key  to  model  residual  useful  life  of  structures  that  are  managed  under  the 
damage  tolerant  approach. 

The  question  is  whether  or  not  these  ideas  would  apply  to  engineering  materials  like  A1 
2024-T351  and  Ti-6A1-4V.  In  the  case  of  the  A1  alloys,  there  is  also  the  issue  of  how  the 
elongated  microstructure  in  the  material  would  affect  the  strain  fields.  Measurements  were 
undertaken  at  several  values  of  AK  at  R=0. 1  How  crack  propagation  direction  with  respect  to  the 
rolling  axis  affected  fatigue  crack  tip  fields  was  investigated,  as  shown  in  Figure  2. 
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Figure  2.  Effect  of  elongated  microstructure  on  fatigue  crack  tip  fields  for  A1  2024-T351.  (a) 
Crack  growth  perpendicular  to  rolling  direction,  (b)  Crack  growth  parallel  to  rolling  direction. 
CT  specimens  tested  under  constant  AK=17  MPa.m0.5  and  R=0.1 


Results  indicated  that  despite  the  clear  difference  in  the  distribution  of  the  strain  fields  for 
the  same  value  of  the  applied  load  (i.e.  the  region  of  large  strain  extends  further  directly  ahead 
of  the  tip  for  case  b)  the  average  crack  growth  rates  were  the  same  within  experimental  error. 
Experimental  results  were  limited  and  full  correlations  between  the  integrated  strain  and  the 
crack  growth  kinetics  could  not  be  obtained,  but  other  information  regarding  crystallographic 
effects  on  fatigue  crack  growth  in  this  material,  e.g.,  the  orientation  of  grains  that  favor  crack 
kinks,  was  obtained.  Nonetheless,  the  results  indicate  that  the  details  of  the  strain  distribution  are 
not  dependent  on  crack  kinetics,  which  favors  non-local  parameters  to  describe  damage. 

Multiscale  damage  criterion 

The  characterized  material  properties  and  observed  crack  nucleation  and  propagation  are 
used  for  the  development  of  multiscale  modeling  and  damage  criteria  in  this  project.  Luo  et  al. 
developed  an  energy  based  damage  evolution  rule  based  on  the  fact  that  fatigue  cracks  in  metals 
tend  to  nucleate  and  propagate  initially  along  slip  planes.  In  this  work,  single  crystal  plasticity  is 
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incorporated  to  the  damage  evolution  rule  so  that  instead  of  calculating  accumulated  fatigue 
damage  along  all  directions  in  3D  space,  calculation  was  carried  out  only  among  potential  active 
slip  systems.  The  corresponding  criteria  are  rewritten  as: 

cr  ( a ) 

dD(“)  =<H2L-  l  >m  (l  +  — — W(Q0  (2) 

CT0  Of 


dY™  =  8onM(d£P)W  + 


1-8 

2 


(js(a)(dyP)(a) 


(3) 


where  m  and  5  are  material  constants.  It  should  be  noted  that  the  direction  of  the  maximum 
material  plane  is  a  unit  vector  along  the  critical  plane  rather  than  a  normal  vector  to  the  plane. 
This  will  be  further  clarified  in  the  results  section.  In  this  work,  all  the  simulations  are  focused 
on  face  centered  cubic  (FCC)  crystal  structures,  such  as  copper  and  aluminum.  A  single  FCC 
crystal  has  12  slip  systems,  comprised  of  four  slip  planes,  each  with  three  slip  directions.  The 
damage  parameter  evolution  of  a  single  grain  of  aluminum  inside  the  meso  RVE  for  12  potential 
active  slip  systems  under  uniaxial  cyclic  loading  was  presented  in  our  previous  work.  The  grain 
was  oriented  for  single  slip  and  the  results  showed  that  the  critical  slip  plane  is  approximately  45 
from  the  loading  direction. 

Using  Eqs.  2  and  3,  the  damage  parameter  D(a)  of  each  slip  system  can  be  calculated.  In 
order  to  incorporate  damage  information  from  the  12  slip  systems,  a  damage  tensor  is  developed 
to  indicate  the  damage  status  at  a  particular  point.  In  Eq.  2,  the  damage  parameters  in  all  12  slip 
systems  are  in  a  strain  energy  density  form.  This  makes  it  reasonable  to  assume  that  the  damage 
tensor,  which  is  also  associated  with  strain  energy  density,  is  a  symmetric  tensor.  On  the  other 
hand,  it  should  reflect  the  directional  effects  for  different  slip  systems.  Thus,  the  relation  between 
the  damage  tensor  increment  and  the  damage  parameter  increment  in  each  slip  system  is 
developed  as  follows: 


dDk  =  nkdDs\  k=l  to  12  (4) 

A  damage  tensor  capable  of  predicting  damage  growth  rate  and  direction  of  damage 
evolution  simultaneously  is  developed  by  computing  the  maximum  eigenvalue  and  the 
corresponding  eigenvector,  and  shown  as  follows: 
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A  simulation  for  a  simple  test  case 
is  illustrated  in  Figure  3.  The  results  are 
slightly  counterintuitive.  At  some  points, 
the  maximum  eigenvalue  becomes  negative, 
which  has  no  physical  meaning. 

Consequently,  the  damage  tensor  is 

decomposed  into  a  deviatoric  part  and  a 
hydrostatic  part.  Since  damage  is  driven  by 
plasticity  and  the  damage  tensor  is 

developed  to  reflect  the  slip  effect,  the 
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Figure  3.  Simulation  for  simple  tension  on  a 


pre-cracked  single  crystal  plate,  (a)  Mesh  of  the 
plate;  (b)  Mises  stress  distribution  in  the  plate; 
(c)  Enlarged  area  at  the  crack  tip. 
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deviatoric  part  of  the  damage  tensor  should  also  contain  the  necessary  information.  The 
maximum  eigenvalue  and  the  corresponding  eigenvector  of  the  deviatoric  part  appear  to  provide 
reasonable  results  in  the  context  of  damage  accumulation. 

Mesoscale  Damage  Vector 


Damage  is  believed  driven  by  plastic  deformation.  Thus,  the  coupling  between  damage 
and  elasticity  is  not  considered.  The  principal  idea  of  this  approach  is  to  develop  a  multiscale 
concept  and  formulate  a  damage  index  for  an  RYE  based  on  current  stress-strain  distribution, 
which  can  represent  the  damage  status  of  the  RVE  considering  some  of  the  microstructure 
features.  Therefore,  several  optimization  methods  and  averaging  techniques  are  used  to  select  the 
critical  local  damage  information  and  transfer  it  to  a  global  damage  variable.  For  this  purpose,  a 
Kreisselmeier-Steinhauser  (KS)  function  based  approach  is  used  to  account  for  the  contribution 
from  all  grains  to  the  total  damage  accumulation  at  mesoscale.  The  KS  function-based  approach 
makes  the  current  multiscale  model  a  statistical  model  rather  than  a  progressive  damage  model. 
It  incorporates  the  criteria  that  allow  contribution  from  the  more  critical  grains  to  be  reflected  in 
the  damage  calculation.  Traditionally,  KS  function  is  used  in  optimization  applications  involving 
multiple  objective  functions  and/or  constraints.  From  a  mathematical  point  of  view,  the  KS 
function  represents  an  envelope  function  for  a  set  of  functions,  as  shown  in  Figure  4  and  defined 
as: 


KS[gi(x)]  =  --ln[Y  e 

P  4^ 

l 

where  p  is  a  parameter  that  determines  the 
proximity  of  the  KS  function  to  the  boundary  of  the 
multiple  objective  functions  h^x).  In  this  work,  the 
multiple  objective  functions  ( ht(x ))  are  the  damage 
growths  of  all  the  grains  in  a  meso  RYE  as  functions 
of  the  time.  When  p  is  positive,  the  KS  function  is 
close  to  the  lower  bound  of  ht(x),  and  when  p  is 
negative,  the  KS  function  is  close  to  the  upper 
bound  of  hi(x).  However,  Eqn.  6  cannot  be  used 
directly  in  this  application  due  to  the  nature  of  the 
exponential  term.  The  modified  KS  function  form  is 
derived  as: 


-phi  (x)i 


(6) 


Figure  4.  KS  function. 


KS[/ifOO]  =  < 


Max(/ij(x))  +  —  ln[^T 
Min  (/if  (x))  — -ln[^ 


e-p(hi(x)-Min  (/ijOO)]! 


upper  bound 
lower  bound 


(7) 


where  the  max  and  min  terms  represent  the  gains  with  the  highest  and  the  lowest  damage 
parameter,  respectively. 


Damage  parameter  and  direction  at  the  mesoscale 

Understanding  the  key  damage  parameters  and  the  crack  direction  are  critical  for  the  robust 
fatigue  modeling  of  metallic  structures.  At  the  mesoscale,  a  1  mmxl  mm  square  is  chosen  as  a 
meso  RVE  (Figure  5)  based  on  Ritchie’s  length  definition  of  a  physically  small  crack  [3].  The 
procedure  to  create  meso  RVE  mesh  comprises  two  steps.  First,  an  Electron  Backscattering 
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Diffraction  (EBSD)  scan  is  used  to  acquire  the 
microstructure  of  the  material  including  grain 
orientations,  grain  shapes,  and  sizes.  Second,  a 
software  package  OOF  (Object-Oriented  Finite 
element  analysis  from  NIST)  is  used  to 
graphically  assign  the  material  properties  to  a 
microstructure  image  for  meshing.  It  should  be 
noted  that  all  the  meso  RVEs  are  generated 
directly  from  the  EBSD  scans  of  the  material  so 
that  all  the  grain  information  used  in  the  FE 
simulation  is  maintained  similar  to  a  real 
microstructure.  Grain  size  and  shape  can  affect 
the  stress-strain  distribution  in  the  RVE,  and  5'  Meso  RVE  containinS  different 

consequently  can  impact  damage  prediction  to  oriented  grains. 

some  extent.  However,  grain  size  and  shape  effects  are  not  explicitly  considered  in  the 
constitutive  model.  Traditional  single  crystal  plasticity  is  used  to  describe  the  material  behavior 
at  the  hotspot  area  of  the  structural  components.  Stress/strain  gradient  effect,  which  accounts  for 
the  size  effect  in  the  constitutive  model  such  as  the  mechanism  based  strain  gradient  crystal 
plasticity  (MSG-CP)  or  other  strain  gradient  theories,  is  not  considered  here.  In  order  to  obtain 
the  damage  parameter  and  direction  at  the  mesoscale,  each  grain  in  the  meso  RVE  (Figure  5)  is 
treated  as  a  single  unit.  The  output  of  each  grain  is  a  damage  vector  D  obtained  by  using  an 
averaging  technique.  Three  steps  are  carried  out  to  calculate  the  damage  vector  in  each  grain. 
First,  the  deviatoric  part  of  the  damage  tensor  is  calculated  for  each  element  within  a  grain.  Next, 
the  deviatoric  damage  tensors  of  all  the  elements  in  the  grain  are  averaged  to  get  a  single 
averaged  damage  tensor  in  which  each  component  of  the  damage  tensor  comes  from  the  mean 
value  of  the  corresponding  components  of  all  the  elements’  deviatoric  damage  tensor.  Finally, 
the  damage  vector  of  each  grain  is  calculated  where  the  magnitude  of  the  damage  vector,  D,  is 
set  equal  to  the  maximum  eigenvalue  of  the  averaged  damage  tensor.  The  direction  of  the 
damage  vector  is  defined  by  the  corresponding  eigenvector. 

The  magnitude  of  the  damage  vector,  D,  in  each  grain  is  a  function  of  time.  The  upper  and 
lower  bound  for  D  of  all  the  grains  within  a  meso  RVE  can  be  obtained  by  applying  the  KS 
function.  The  damage  index  for  the  meso  RVE  can  be  defined  as: 

Dmeso  =  0 9KSU  +  (1  -  6)KSi)/Dc  (8) 

where  0  is  related  to  the  critical  damage  value  and  the  total  number  of  grains  within  a  meso  RVE. 
Physically,  this  parameter  measures  the  number  of  grains  that  reaches  the  critical  damage  value, 
Dc  .  The  critical  damage  value  is  also  used  to  determine  crack  initiation  in  the  meso  RVE.  The 
damage  direction  in  the  meso  RVE  can  be  obtained  by  normalizing  the  sum  of  all  damage 
vectors  for  all  grains.  For  prediction  of  fatigue  crack  initiation,  the  criteria  should  be  related  to 
the  local  damage  parameter,  i.e.,  the  damage  parameter  D  of  each  grain.  When  the  maximum 
damage  parameter  of  each  grain  within  a  meso  RVE  reaches  the  critical  damage  threshold,  the 
corresponding  number  of  fatigue  cycles  is  treated  as  fatigue  crack  initiation  and  the 
corresponding  grain  with  the  maximum  damage  parameter  is  regarded  as  the  crack  initiation 
location.  Moreover,  the  goal  of  using  a  damage  index  for  the  meso  RVE,  Dmeso,  is  to  determine 
the  failure  of  the  RVE  by  checking  whether  Dmeso  reaches  a  value  of  one  or  not.  In  this  work,  for 
convenience  of  experimental  validation,  the  failure  crack  length  of  the  meso  RVE  is  taken  as 
1mm  considering  the  size  of  the  meso  RVE  and  the  resolution  of  the  digital  image  acquisition 
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system  used  in  the  experiments.  The  meso  RVE  size  is  about  lmmxlmm  and  only  the  failure  of 
the  meso  RVE  is  validated  considering  our  current  experimental  capabilities. 

Development  of  statistical  volume  element  model 

In  order  to  efficiently  implement  the  multiscale  analysis  and  prediction  of  fatigue  life,  it 
is  necessary  to  design  a  model  to  represent  material  properties  and  mechanical  responses  for 
polycrystalline  materials.  The  primary  goal  of  developing  an  SVE  model  is  to  improve  the 
computational  efficiency  and  reduce  the  associated  preprocessing  effort  while  maintaining 
accuracy.  This  is  achieved  by  simplifying  grain  shapes 
in  the  SVE  model,  which  results  in  ease  of  assembly, 
reduction  in  preprocessing  time,  and  reduction  in  total 
number  of  elements  (irregular  grain  shapes  require  more 
and  finer  elements).  The  SVE  model  is  constructed  by 
assembling  grains  whose  features  are  statistically 
sampled  from  a  pool  of  measured  experimental 
characterization  data  using  Electron  Backscatter 
Diffraction  (EBSD)  scan  methods  (Figure  6a).  This 
approach  provides  a  computationally  efficient 
alternative  to  traditional  techniques  while  maintaining 
simulation  accuracy. 

The  key  features  of  grains  include  grain 
orientation,  misorientation,  size,  shape,  aspect  ratio,  and 
principal  axis  direction.  Previous  statistical  analysis  has 
proven  that  there  is  no  obvious  trend  in  distribution  of 
grain  orientations  in  aluminum  alloy  2024-T351. 

Therefore,  grain  orientations  are  assumed  to  be 
randomly  distributed.  The  grain  orientations  used  in  the 
SVE  model  are  randomly  selected  from  an  experimental 
dataset.  The  distribution  of  misorientation  is  chosen  to 
closely  match  that  of  actual  data  from  experimental 
observation.  For  A1  2024-T351,  the  grain  size  in  the 
SVE  model  is  0.0222  mm2.  It  is  close  to  the  average 
measured  value  from  the  observation,  which  is  0.0225 
mm2.  This  preserves  the  actual  density  of  grains  in  the  SVE  model.  A  rectangular  grain  shape  is 
selected  to  easily  assemble  the  grains  within  the  square  SVE.  The  aspect  ratio  is  set  to  be  0.4, 
which  is  close  to  the  actual  average  measured  value  of  0.384.  The  principal  axis  direction  is 
chosen  to  be  the  same  as  the  rolling  direction  of  the  material.  The  final  architecture  of  the  SVE  is 
shown  in  Figure  6  b).  Figure  7  shows  three  pole  figures  for  grain  orientations.  The  distributions 
do  not  exhibit  any  trends,  and  therefore  are  assumed  to  be  normally  distributed.  It  is  noted  that 
when  formulating  the  SVE,  the  orientations  were  chosen  at  random  from  the  pool  of  available 
material  characterization  data.  The  misorientation  is  determined  from  the  arrangement  of  the 
grains  in  the  SVE.  Since  the  ordering  can  be  controlled,  the  misorientation  was  chosen  to  closely 
match  that  of  the  experimental  data  shown  in  Figure  8  a).  To  accomplish  this,  nine  randomly 
orientated  SVEs  were  generated.  A  typical  histogram  of  their  misorientation  is  shown  in  Figure  8 
b). 


Figure  6.  Comparison  between:  a) 
actual  microstructure  scan;  b) 
constructed  SVE.  Colors  represent 
different  grain  orientations. 
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Figure  7  .  Flattened  3D  polar  plots  of  crystal  orientation:  a)  {1,0,0};  b)  { 1,1,0};  c)  { 1,1,1 }. 


Figure  8.  Statistical  distribution  used  for  misorientation  feature  identification  a)  actual  data;  b) 
example  of  generated  SVE. 

To  validate  the  SVE  based  multiscale  modeling  approach,  two  measures  were  considered: 
hot  spot  location  and  far  field  distributions.  To  ensure  the  SVE  is  representative  of  the  bulk 
material  and  mechanical  responses,  both  measures  must  remain  consistent  with  or  without  an 
SVE  implemented  into  the  finite  element  model.  The  SVE  is  expected  to  have  a  direct  influence 
on  the  local  fields  near  the  hot  spot  due  to  local  crystal  plasticity  effects.  Without  this,  there 
would  be  no  means  to  calculate  the  damage  at  the  microscale  for  fatigue  life  estimation.  Figure  9 
shows  a  comparison  of  von  Mises  stress  distributions  in  the  lug  joint  models,  with  and  without 
an  SVE.  The  overall  contours  show  similar  trends,  i.e.,  the  implemented  SVE  has  a  small  effect 
on  the 

mechanical 
response.  The 
hot  spot 
location 
remains 
unchanged  at 
the  inside 
corner  radius 
below  the  a) 

loading  pins  Figure  9.  Local  von  Mises  stress  field  distribution  within  Lug  Joint  a)  with  SVE; 
and  the  b)  without  SVE. 


magnitude  remains  approximately  the  same.  The  distribution  of  stress  in  the  far  field  outside  of 
the  hot  spot  position  also  remains  the  same.  With  both  the  hot  spot  and  far  field  values  nearly 
unchanged,  the  SVE  based  multiscale  modeling  technique  is  proven  capable  of  representing  the 
homogeneous  material. 

Simulation  and  experiment  results:  fatigue  tests  on  lug  joints 

Uniaxial  fatigue  tests  using  lug  joint  specimens  are  conducted  to  validate  predicted 
fatigue  crack  growth  estimated  by  the  developed  SVE  model  and  the  damage  criterion.  Lug  joint 
specimens  are  machined  from  a  bulk  A1  2024-T35 1  plate.  The  rolling  direction  is  the  same  as  the 
lengthwise  direction  of  the  lug  joint  specimen.  The  specimens  are  cyclically  loaded  in  a  servo 
hydraulic  desktop  test  frame  using  a  sinusoidal  waveform  with  a  load  range  between  489.3  N  and 
4893.0  N  (load  ratio  of  0.1)  at  a  rate  of  20 
Hz.  To  validate  the  multiscale  model,  the 
number  of  cycles  necessary  to  achieve  a  1 
mm  crack  is  obtained  experimentally.  A  total 
of  five  specimens  are  tested.  The 
experimental  results  are  summarized  in  Table 
1.  Figure  10  shows  the  tested  specimen  with 
an  illustrated  crack  angle  in  the  defined 
coordinate  system.  The  number  of  cycles  to 
produce  a  1  mm  crack  ranges  from  68  K  to 
98  K  cycles,  and  the  crack  direction  is  around 
-4°  or  -28°. 

Table  1.  Experimental  results  of  fatigue  life  and  crack  direction. 


Samplel  Sample2  Sample3 


Sample4  Sample5 


Figure  10.  Crack  formation  in  experimental 
samples. 


Sample 

1 

2 

3 

4 

5 

Average 

Fatigue  life  (No.  of  cycles) 

92K 

98K 

68K 

83K 

82K 

84.6K 

Crack  direction  (degrees) 

-25 

-34 

-29 

-28 

-28 

-28.8 

For  computational  efficiency,  a  two-scale  mesh  is  used  in  the  multiscale  analysis  of 
complex  structural  components  such  as  a  lug  joint.  Preliminary  stress  analysis  has  been 
conducted  by  applying  homogeneous  elasto-plastic  material  model  to  identify  the  hotspot  of  the 
structural  component  in  ABAQUS.  The  meso  RVE  mesh  generated  using  the  software  OOF  has 
been  used  at  the  hotspot  of  the  lug  joint.  The  rest  of  the  lug  joint  is  described  as  a  homogenous 
material.  The  constitutive  relation  of  this  homogeneous  material  was  obtained  by  homogenizing 
the  meso  RVE  stress-strain  response.  First,  a 
force  was  applied  at  the  right  edge  of  the  meso 
RVE  where  plane  stress  elements  were  used. 

The  UMAT  based  on  single  crystal  plasticity 
has  been  used  to  describe  material  behavior  for 
each  grain  within  the  meso  RVE.  Then,  the 
displacement  of  each  node  at  the  edge  was 
calculated.  The  corresponding  homogenized 
meso  RVE  stress-strain  response  can  be  plotted 
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based  on  the  applied  load  and  the 
displacement  of  each  node.  The  two- 
scale  mesh  was  generated  using  the 
commercial  software  Altair  Hypermesh. 

Figure  1 1  shows  the  two-scale  mesh  of 
the  lug  joint.  For  this  work,  all  the 
simulations  were  carried  out  in  2D. 

Three  nodes  and  four  nodes  plane 
stress  elements  (CPS3  and  CPS4)  are 
used  for  the  FE  simulations  of  the  lug 
joint  to  investigate  the  surface  of  the 
specimen.  Symmetric  boundary 
conditions  were  used  for  simulation  so 
that  only  half  of  the  lug  joint  was  analyzed  in  ABAQUS.  The  same  cyclic  loading  condition 
which  was  used  during  the  fatigue  test,  that  is,  490  N  (110  lbs)  to  4900  N  (1100  lbs)  with  a 
frequency  of  20  Hz  in  sinusoidal  waveform  applied  at  the  pin  hole  of  the  lug  is  used  in 
simulation.  Figure  12  shows  the  von  Mises  stress  distribution  in  the  lug  joint  under  simple 
tension  with  the  enlarged  hotspot  area  showing  the  non-uniform  distribution  due  to  different 
grain  orientations. 

The  developed  damage  tensor  was  implemented  in  the  UMAT  and  a  data  processing  code 
was  also  developed  in  Matlab.  Damage  evolution  for  all  grains  in  the  meso  RVE  is  plotted  in 
Figure  13  for  20  cycles.  Figure  13  shows  that  after  10  cycles,  the  damage  evolution  in  each  grain 
becomes  stable,  which  provides  a  basis  for  using  a  linear  fit  to  extrapolate  the  damage  evolution 
in  individual  grains.  It  should  be  noted  that,  the  focus  of  this  work  was  to  propose  a  new 
methodology  for  fatigue  damage  prediction  taking  into  consideration  microstructure  features. 
Therefore,  all  simulations  conducted  and  presented  here  are  under  constant  cyclic  loading.  Thus, 
the  damage  evolution  in  each  grain  becomes  almost  linear  after  10  cycles.  For  random  loading 
conditions,  future  work  will  address  building  a  relationship  between  the  applied  load  and  the 
damage  growth  in  individual  grains.  Figure  14  a)  shows  the  damage  parameter  of  all  the  grains 
within  the  RVE  at  time,  t=l. 2667s,  and  Figure  14  b)  presents  the  enlarged  area  of  the  RVE  and 
highlights  the  crack  initiation  area  where  the  damage  parameter  was  maximum.  This  critical 
grain  is  labeled  as  grain  no.  9,  shown  in  Figure  14  a).  The  results  confirm  that  the  critical  grain  in 
the  RVE  is  located  close  to  the  free  surface  at  the  shoulder  of  the  lug  joint. 


Figure  12.  Mises  stress  distribution  of  lug  joint. 
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a) 

Figure  13.  a)  Damage  evolution  in  each  grain  for  20  cycles,  b)  An  example  shows  that  damage 
grows  linearly  after  10  cycles  in  each  individual  grain. 
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A  Matlab  program  was  developed  for 
post-processing  of  finite  element  results. 
Figure  15  shows  the  microstructure  and  the 
grain  size  distribution  of  the  meso  RVE  used. 
Figure  16  presents  the  normalized  damage 
parameter  for  the  meso  RVE  versus  the 
number  of  cycles.  The  simulation  data  fit  a 
quadratic  polynomial  well,  which  is  intuitive 
given  that  as  more  grains  reach  the  critical 
damage  value,  the  accumulation  of  damage  in 
the  meso  RVE  accelerates.  The  estimated 
number  of  cycles  until  failure  in  the  meso 
RVE,  i.e.,  when  the  normalized  damage  index 
of  the  meso  RVE  reaches  an  unit  value,  is 
208Kcycles.  The  result  of  the  simulation 
matches  well  with  the  experimental  results  of 
samples  1-3  shown  in  Table  2.  The 
corresponding  eigenvector,  an  indicator  of  the 
potential  damage  direction,  is  calculated  by 
the  weighted  average  method  from  all  the 
grains.  A  histogram  showing  the  frequency  of 
damage  occurrence  along  a  particular 
direction  is  presented  in  Figure  17  a).  The 
figure  indicates  that  the  directions  of 
maximum  damage  in  the  RVE,  obtained  from 
simulation,  are  approximately  -30°and  52  °. 
The  experimental  crack  directions  from  lug 
joint  fatigue  tests  are  shown  in  Figure  17  b). 
Comparing  the  simulation  results  with  the 
experimental  data,  one  of  the  potential  damage 
directions  (-30°)  obtained  from  the  simulation 
matches  the  experiments  well.  Further 
simulations  were  conducted  to  verify  the 
model  and  consider  uncertainty  of  the  meso 
RVE.  The  same  lug  joint  with  only  one  meso 
RVE  located  at  the  shoulder  was  used.  The 
meso  RVEs,  however,  contain  different 
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Figure  14.  a)  Damage  distribution  in  RVE  for 
all  grains,  b)  Location  of  critical  grain. 


Figure  15.  Microstructure  and  grain  size 
distribution  of  the  meso  RVE 


Figure  16.  Normalized  damage  index  for 
meso  RVE  vs.  No.  of  cycles 


oriented  grains  and  different  number  of  grains. 

As  mentioned  before,  all  meso  RVE  sizes  are  approximated  as  lmmxlmm.  The  average  grain 
size  and  estimated  failure  of  meso  RVEs  are  shown  in  Table  2.  It  can  be  observed  that  the 
estimation  of  RVE  failure  in  simulation  5  is  close  to  the  fatigue  test  results  obtained  in  sample  4. 
Variability  in  the  simulation  results  indicates  that  the  fewer  number  of  cycles  required  in  sample 
4  to  obtain  a  1mm  crack  could  be  due  to  different  oriented  grains. 

The  five  simulations  suggested  that  grain  orientation  will  affect  the  estimation  of  RVE 
failure.  The  common  feature  in  all  five  simulations  is  that  the  RVEs  used  in  the  simulations  are 
generated  from  the  scans  taken  from  the  same  A1  alloy  plate.  However,  those  scans  are  not 


12 


directly  taken  from  the  lug  joint  samples,  which 
results  in  the  variability  of  the  input  data  for  the 
model.  To  prevent  this  variability,  an  RVE  which 
is  directly  scanned  from  the  lug  joint  shoulder  is 
created.  The  procedure  involves  obtaining  four 
scans  from  both  shoulders  of  the  lug  joint  on  both 
sides  and  conducting  a  fatigue  test  on  the  same 
lug  joint  sample  to  determine  which  scan  should 
be  used  in  the  validation  simulation.  The  lowest 
number  of  cycles  recorded  was  33k  cycles  to 
initiate  a  1.5  mm  crack. 

The  simulation  result  of  RVE  failure  is 
found  to  be  12.4  k  cycles.  Considering  the 
number  of  cycles  required  for  short  crack 
propagation  up  to  1.5  mm  that  can  be  obtained 
experimentally,  the  author  believes  the  simulation 
result  is  acceptable.  The  direction  of  crack 
propagation  at  the  early  stage  of  fatigue  test  is 
around  -53°  with  respect  to  the  horizontal 
direction.  Results  show  two  major  potential  crack  direction  bands.  One  is  from  approximately  - 
58°  to  -43°  and  the  other  is  from  60°  to  85°.  The  experiment  result  shows  the  crack  propagated 
along  one  of  the  potential  crack  direction  ranges  obtained  from  the  simulation,  indicating  that  the 
model  can  predict  potential  crack  directions. 

Table  2.  Average  grain  size  and  estimation  of  failure  for  different  meso  RVEs 


Lug  Joint 
Simulations 

Simulation 

1 

Simulation 

2 

Simulation 

3 

Simulation 

4 

Simulation 

5 

Average  grain  size 
(pm) 

208.591 

191.866 

190.266 

242.684 

191.082 

Estimate  failure  of 
RVE  (K  cycles) 

208 

184 

298 

177 

116 
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Figure  17.  a).  Histogram  of  damage 

direction  in  RVE;  b).  Cracking  directions 
from  fatigue  tests. 


Task  2:  Methods  for  In- Situ  Interrogation  and  Detection 

Damage  detection  and  characterization  is  critical  for  the  development  of  an  integrated 
SHM  system.  An  adaptive  learning  baseline-free  structural  damage  estimation  method  has  been 
developed  that  is  robust  to  variability  in  environmental  and  operating  conditions.  Physics-based 
damage  growth  models  have  been  integrated  with  sensor  signal  processing  algorithms  for 
estimating  progressive  structural  damage  effectively  using  sequential  Monte  Carlo  techniques. 
Within  this  framework,  performance  bounds  have  been  evaluated  and  optimal  sensor  scheduling 
and  automatic  measurement  screening  approaches  were  designed  and  implemented  to  further 
maximize  the  damage  estimation  performance.  The  applicability  of  compressive  sensing 
techniques  was  studied  for  substantially  reducing  structural  data  volume  in  order  to  save  on 
storage,  transmission,  and  processing  costs.  A  transfer  learning  methodology  was  proposed  that 
utilizes  knowledge  transfer  between  multiple  distributed  sensors  in  order  to  increase  structural 
damage  classification  accuracy.  Information  obtained  from  a  deterministic  physics  based 
multiscale  model  was  utilized  to  formulate  a  new  statistical  model  for  describing  small  crack 
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growth  in  aluminum  2024  T351.  Vibration  based  and  laser  ultrasonic  based  methods  were  also 
studied  in  this  project. 

Adaptive  learning  based  structural  damage  estimation 

A  novel  adaptive  learning  structural  damage  detection  and  estimation  methodology  has 
been  developed  for  robustness  to  variable  or  changing  conditions.  Variable  environmental  and 
operating  conditions  are  known  to  be  a  major  challenge  in  the  design  of  reliable  structural  health 
monitoring  (SHM)  systems  because  both  damage  evolution  and  sensor  data  are  often  strongly 
influenced,  for  example,  by  changes  in  temperature,  geometry  or  configuration.  Unlike 
conventional  statistical  methods  that  can  be  limiting  in  highly  complex  and  rapidly  changing 
environments,  the  adaptive  learning  approach  uses  advanced  stochastic  modeling  to  continuously 
evolve  with  the  time-varying  environment. 

Specifically,  Dirichlet  process  (DP)  mixture  models  are  utilized  to  automatically  adapt  to 
structure  within  the  data.  The  DP  mixture  model  provides  for  a  growing,  possibly  infinite 
number  of  mixture  components  or  clusters/classes,  a  finite  number  of  which  are  manifested 
within  the  given  data  (time-frequency  probability  density  function  (PDF)  statistical  distance 
features  extracted  from  sensor  signals).  The  appropriate  number  of  mixture  components  and 
mixture  proportions,  and  the  mixture  distribution  parameters,  are  learned  adaptively  from  the 
data.  The  learning  of  DP  mixture  model  parameters  is  performed  efficiently  using  Markov  Chain 
Monte  Carlo  (MCMC)  techniques.  The  adaptively  identified  components  or  classes  can  then  be 
linked  to  different  types  of  damage  within  a  structure  or  different  possible  variations  in  the 
external  environment  for  the  same  type  of  damage.  The  damage  state  inference  is  performed 
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Figure  18.  Block  diagram  of  the  adaptive  learning  based  structural  damage  estimation 
method. 


using  a  Bayesian  filter  that  combines  the  adaptive  data  model  with  a  physics  based  progressive 
damage  model.  The  main  advantage  of  this  approach  is  that  no  baseline  training  data  is  required 
and  signals  can  be  classified  on  the  fly  to  new  (previously  unseen)  damage  classes,  yielding  an 
adaptive  and  effective  approach  for  online  SHM.  Figure  18  shows  a  block  diagram  of  the 
adaptive  learning  based  structural  damage  estimation  method. 

The  proposed  adaptive  learning  technique  has  been  validated  for  the  detection  of  fatigue 
damage  in  metallic  structures  under  variable  conditions.  Figure  19  a)  shows  a  comparison  of  the 
true  and  estimated  crack  lengths  using  Lamb  wave  measurements  from  an  aluminum  compact 
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tension  (CT)  specimen  subjected  to  variable-amplitude  fatigue  loading.  Figure  19  b)  shows  the 
true  and  estimated  fatigue  crack  lengths  using  impedance  measurements  from  an  aluminum  lug 
joint  specimen  under  varying  temperature  conditions.  The  results  show  good  damage  estimation 
performance  and  robustness  to  variability.  It  should  be  noted  that  in  the  lug  joint  experiments, 
the  temperature  was  changed  only  during  sensor  data  acquisition  and  not  while  actual  fatiguing 
of  the  sample.  More  extensive  experiments  are  underway  to  test  the  adaptive  method’s 
performance  under  realistic  temperature  variability  scenarios. 
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Figure  19.  Adaptive  learning  results  for  fatigue  crack  damage  estimation  in  metallic  structures 
under  variable  conditions,  (a)  aluminum  CT  specimen  subjected  to  variable-amplitude  loading, 
(b)  aluminum  lug  joint  specimen  under  varying  temperature  conditions. 


Progressive  structural  damage  estimation  using  sequential  Monte  Carlo  techniques 

A  new  state-space  framework  based  method  has  been  developed  for  progressive 
structural  damage  estimation  using  sequential  Monte  Carlo  techniques.  The  method  integrates 
physics-based  models  with  information  from  sensor  measurements  for  effectively  estimating  and 
tracking  an  evolving  damage  state  using  sensor  data  collected  online.  The  state-space 
formulation  is  used  to  model  (a)  the  damage  evolution  process  (assumed  Markov),  and  (b)  the 
measurement-damage  relationship.  The  damage  state  estimation  is  then  carried  out  in  a  Bayesian 
setting.  When  the  state  and  measurement  dynamics  are  linear  and  the  noise  can  be  modeled  as 
Gaussian,  the  solution  can  be  computed  using  the  Kalman  filtering  technique.  However,  in  the 
damage  estimation  problem  the  underlying  physical  process  is  generally  nonlinear  and  non- 
Gaussian,  and  the  tracking  problem  is  solved  using  sequential  Monte  Carlo  techniques,  e.g., 
particle  filter  (PF).  In  the  particle  filter,  the  desired  posterior  probability  distribution  over  the 
state  variables  is  summarized  using  samples  (particles  and  weights).  This  representation  is  then 
updated  continuously  to  include  the  information  provided  by  recent  measurements. 

The  proposed  method  was  implemented  for  the  estimation  of  fatigue-induced  crack 
damage  in  an  aluminum  CT  specimen  under  cyclic  loading.  The  system  model  was  obtained 
from  fracture  mechanics  based  on  Paris’  Law.  The  stochastic  relationship  in  the  measurement 
model  was  obtained  from  hidden  Markov  modeling  (HMM)  of  joint  time-frequency  (TF) 
features  extracted  from  the  PZT  sensor  signals  using  the  Matching  Pursuit  Decomposition  (MPD) 
algorithm.  Figure  20  shows  results  of  estimating  the  crack  length  in  the  CT  specimen  using  the 
proposed  particle  filtering  method  with  discrete  and  continuous  measurement  density  HMMs.  It 
can  be  seen  that  the  algorithm’s  performance  is  very  good,  with  only  a  few  millimeters  of  error 
in  the  estimated  crack  lengths  with  respect  to  the  experimentally  observed  values. 
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^  N  unber  of  fatigue  cycles  (Idb-cyc)  ^  Nunberof  fatigue  cycles  (tib-cyc) 

Figure  20.  Fatigue  crack  estimation  in  a  CT  specimen  using  sequential  Monte  Carlo  (particle 
filtering),  (a)  discrete  HMMs,  (b)  continuous  HMMs. 

Compressive  sensing  of  structural  data 

Compressive  sensing  (CS)  techniques  have  been  investigated  in  order  to  reduce  the 
amount  of  structural  monitoring  data  collected  for  achieving  savings  in  storage,  transmission  and 
processing  related  costs.  Specifically,  it  was  observed  that  the  sparsity  of  the  dispersive  Lamb 
wave  measurements  in  time-frequency  and  time-scale  bases  can  be  exploited  to  substantially 
reduce  data  volume  at  the  sensing  stage.  In  our  experiments,  we  considered  PZT  sensor  signals 
recorded  from  a  bolted  aluminum  plate.  About  only  20%  of  the  samples  of  a  signal  were  first 
randomly  selected  as  measurements  and  these  compressively  sensed  noisy  measurements  were 
then  decoded  by  solving  a  Basis  Pursuit  Denoising  convex  optimization  problem  (with  a  wavelet 
packet  basis  used  for  the  signal  representation).  It  was  shown  that  the  structural  signals  could  be 
accurately  reconstructed  for  subsequent  processing,  as  see  Figure  21  below. 


Figure  21.  Compressive  sampling  of  a  PZT  sensor  signal  recorded  from  a  bolted  aluminum  plate, 
(a)  original  complete  signal  and  20%  randomly  selected  samples,  (b)  original  signal  and  CS 
recovered  approximation. 

Adaptive  measurement  selection  for  progressive  damage  estimation 

Noise  and  interference  in  sensor  measurements  degrade  the  quality  of  data  and  have  a 
negative  impact  on  the  performance  of  structural  damage  diagnosis  systems.  We  have  designed  a 
novel  adaptive  measurement  screening  algorithm  to  automatically  select  the  most  informative 
sensor  measurements  and  use  them  for  fast  and  accurate  structural  damage  estimation  within  the 
SMC  framework.  Our  approach  is  motivated  by  the  data  association  concept  which  has  been  well 
studied  and  used  successfully  in  the  problem  of  radar  based  target  tracking.  The  idea  is  to 
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adaptively  screen  the  sensor  measurements  based  on  their  plausibility  given  the  most  recent 
damage  state  estimates  and  previous  measurement  history.  This  is  determined  by  evaluating  the 
measurement  likelihoods  under  the  predicted  measurement  distribution. 

Two  approaches  were  considered  for  realizing  the  adaptive  measurement  selection 
method:  the  Maximum  Likelihood  (ML)  approach  and  the  Minimum  Mean-Squared  Error 
(MMSE)  approach.  The  ML  approach  selects  the  current  single  most  likely  measurement  from 
the  set  of  all  available  measurements  and  uses  it  to  estimate  the  damage  state.  The  MMSE 
approach  computes  the  damage  state  estimate  as  a  weighted  average  of  the  estimates  obtained 
using  each  of  the  measurements  individually,  with  the  weights  being  proportional  to  the 
likelihood  of  the  respective  measurements;  unlike  the  ML  approach,  it  puts  more  emphasis  on 
measurements  which  are  more  likely  and  allows  information  from  multiple  measurements  to 
contribute  to  the  damage  state  estimate.  Both  approaches  effectively  reject  noisy  and/or 
implausible  measurements,  and  in  doing  so,  screen  for  outliers  and  retain  only  the  most 
informative  measurements. 

The  proposed  approach  was  applied  to  further  increase  the  performance  of  the  particle 
filtering  progressive  fatigue  damage  estimation  method  for  tracking  crack  length  in  an  aluminum 
CT  specimen  using  noisy  PZT  sensor  measurements.  Discrete  HMMs  were  utilized  to  model  TF 
features  extracted  from  the  sensor  signals  using  the  MPD  algorithm.  Figure  22  shows  the  results 
obtained  with  noise  suppression  using  ML  and  MMSE  based  adaptive  measurement  selection 
approaches  and  using  randomly  selected  single  and  multiple  measurements.  It  can  be  seen  that 
suppressing  the  use  of  noisy  measurements  with  the  proposed  MMSE  based  adaptive 
measurement  selection  algorithm  helped  to  significantly  improve  the  damage  estimation 
performance,  demonstrating  up  to  80%  reduction  in  the  MSE  for  estimating  crack  length.  The 
ML  based  adaptive  measurement  selection  employs  hard  decisions  and  is  somewhat  sensitive  to 
prediction  errors;  in  this  example  it  did  not  deliver  the  as  much  improvement  in  the  performance 
as  the  MMSE  approach. 
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Figure  22.  Progressive  damage  estimation  in  a  CT  sample  without  and  with  noise  suppression 
using  adaptive  measurement  selection,  (a)  crack  length  estimation  results,  and  (b)  comparison  of 
mean  squared  error. 


Structural  damage  classification  with  insufficient  data  using  transfer  learning 

In  many  structural  damage  classification  scenarios,  a  large  amount  of  ‘training  data’  may 
be  required  for  learning  statistical  parameters,  which  is  often  difficult  to  obtain  and  leads  to 
compromised  system  performance  under  these  data-scarce  conditions.  However,  in  many 
applications  a  modest  amount  of  data  may  be  available  from  a  few  different  but  related 
experiments.  We  have  developed  a  new  structural  damage  classification  method  that  makes  use 
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Figure  23.  Transfer  Learning  domains  used  for 
fatigue  damage  classification  in  the  lug 
specimen. 


of  statistics  from  related  tasks  to  improve 
the  classification  performance  on  a  data  set 
with  limited  training  examples.  The 
approach  is  based  on  the  framework  of 
transfer  learning  (TL),  which  provides  a 
mechanism  for  knowledge  transfer  between 
related  learning  tasks  (from  a  ‘source’ 
domain  to  a  ‘target’  domain)  and  thereby 
reduces  the  training  burden  without 
significantly  compromising  on  the  learner’s 
performance. 

The  proposed  transfer  learning  based  damage  classification  method  shares  model 
statistics  learned  from  multiple  distributed  sensors  to  achieve  significant  performance  increase  in 
classifying  fatigue  damage  in  an  aluminum  lug  joint  specimen  using  a  limited  amount  of  PZT 
sensor  data.  Specifically, 
data  was  collected  using 
four  sensors,  but 
sufficient  training  data 
was  available  only  from 
one  sensor  which  was 
considered  as  the  source 
domain.  The  remaining 
three  sensors  with 
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(e)  Classification  of  sensor  4  signals  (without  TL). 
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(f)  Classification  of  sensor  4  signals  (with  TL). 


Figure  24.  Fatigue  damage  classification  in  an  aluminum  lug  joint 
specimen  using  PZT  sensor  data  without  and  with  the  use  of 
Transfer  Learning. 


18 


each  class  for  training.  On  the  other  hand,  only  5  signals  were  used  for  training  in  each  of  the 
target  domains.  The  proposed  TL  based  structural  damage  classification  method  was  then 
applied.  Figure  24  shows  a  comparison  of  the  fatigue  damage  classification  results  without  and 
with  the  use  of  TL.  Our  results  confirm  the  penalty  of  using  inadequate  amount  of  training  data 
for  the  damage  classification.  As  seen  from  the  plots,  there  is  a  dramatic  improvement  in  damage 
classification  performance  when  TL  is  used  to  share  knowledge  between  the  source  and  target 
domain  sensors. 

Integration  of  physics  based  multiscale  models  into  progressive  damage  estimation  framework 

The  growth  behavior  of  microstructure  level  cracks  (cracks  of  size  smaller  than  1mm)  is 
different  from  that  of  macro  level  cracks  and  cannot  be  accurately  modeled  using  the  simple 
Paris’  Law  relationship.  While  several  deterministic  small  crack  growth  models  exist  in  the 
literature,  a  probabilistic  model  is  desirable  in  order  to  be  able  to  account  for  the  uncertainty  in 
crack  growth  due  to  randomness  at  the  microstructure  level.  In  our  work,  data  on  crack  growth  in 
the  micro  length  scales  obtained  from  the  physics  based  multiscale  model  developed  in  Dr. 
Chattopadhyay’s  group  is  being  utilized  to  design  a  new  statistical  model  for  describing  small 
crack  growth  in  aluminum  2024  T351. 

The  proposed  statistical  model  characterizes  the  small  crack  growth  process  using  two 
distinct  growth  states:  normal  growth  and  accelerated  growth.  The  normal  and  accelerated 
growth  states  are  modeled  using  linear  growth  regimes  with  independent  truncated  Gaussian  and 
Gamma  distributed  random  variables  for  crack  growth  rate  and  fatigue  life  cycles.  The 
parameters  in  the  model  are  estimated  with  maximum-likelihood  estimation  (MLE)  using  data 
obtained  from  FEM  simulations  of  the  physics  based  multiscale  model  with  several  A1  2024 
T351  microstructure  instantiations.  The  model  can  be  adapted  for  use  within  the  sequential 
Monte  Carlo  based  progressive  damage  estimation  framework  in  the  small  crack  growth  regime. 

Vibration-based  health  monitoring  with  temperature  variations 

Vibration-based  damage  detection  has  been  considered  as  a  useful  tool  for  damage 
detection,  classification,  and  in-situ  monitoring.  Vibration  based  damage  detection  methods 
considering  the  temperature  effects  have  been  studied  in  this  project.  The  beam  studied  here,  as 
shown  in  Figure  25  a),  is  assumed  to  be  a  uniform  Euler-Bemoulli  beam  actuated  with  a  PZT 
actuator,  undergoing  temperature  variations.  It  is  also  assumed  to  have  a  single  fatigue  crack. 
The  crack  is  treated  as  a  localized  reduction  in  the  stiffness  and  modeled  as  a  massless  rotational 
spring  at  the  location  of  the  crack.  The  beam  is  then  considered  to  be  of  bisections  connected  by 
this  spring,  as  shown  in  Figure  25  b).  The  beam  is  assumed  to  have  a  length  L,  width  b  and 
thickness  hb  and  a  uniform  cross  section  along  the  beam.  For  simplicity  in  the  modeling,  the 
single  fatigue  crack  is  assumed  to  be  non-breathing  during  the  deformation  of  the  beam.  The 
surface-bonded  piezoceramic  (PZT)  actuator  is 
assumed  to  have  the  length  of  La ,  width  of  b  and 
thickness  of  ha,  attached  at  the  xla  from  the  start  of 
the  beam. 

Change  in  the  temperature  distribution  inside 
the  beam,  T(x,y,z,t)  ,  affects  the  vibratory 
characteristics  of  this  system  through  alterations  of 
material  properties,  like  modulus  of  elasticity,  as 
well  as  causing  thermal  stresses  inside  the  beam. 

Thermal-induced  changes  of  the  dimensions  of  the 


Figure  25.  Beam  with  a  single  fatigue 
crack  under  PZT  actuation  and 
temperature  variations. 
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beam  through  material  expansionis  negligible  compared  to  the  effects  on  the  material  properties 
and  thermal  stresses  and  therefore  have  not  been  considered  in  the  present  modeling. 

Temperature-dependent  shift  in  the  frequency  is  tabulated  for  an  Euler-Bernoulli  beam  of 
length  9.45",  width  0.75"  and  thickness  0.18",  made  of  aluminum  (mass  density  of  168.56 
lb/ft3(2700  kg/m3))  in  Table  3.  Temperature  is  varied  from  0°C  (32°F)  to  220°C  (428°F)  and  the 
first  natural  frequency  of  the  beam  is  calculated  for  the  changing  modulus  of  elasticity  in  this 
range  of  temperature.  These  frequencies  are  calculated  based  on  simple  vibrations  modeling  of 
an  Euler-Bernoulli  beam,  considering  only  the  first  mode  of  vibrations. 

Table  3.  Shift  in  the  first  natural  frequency  of  an  aluminum  beam  due  to  the  temperature-varying 


modulus  of  elasticity  for  c 

ifferent  boundary  conditions. 

Boundary  Condition 

Cantilever 

Pinned-Pinned 

Clamped-Pinned 

Clamped-Clamped 

Frequency  Shift  (Hz/°C) 

0.0188 

0.0529 

0.0828 

0.1199 

For  a  specific  input  voltage  applied  on  the  PZT  and  a  known  temperature  distribution, 
generalized  coordinates  of  the  system  will  be  derived  once  the  matrices  of  mass,  stiffness,  PT  (t) 
mdrj  are  formulated.  The  values  of  these  matrices  depend  on  the  material  properties  of  the  beam 
(which  are  functions  of  the  temperature  inside  the  beam)  and  the  thermal  moments  inside  the 
beam,  as  well  as  the  choice  of  admissible  functions  of  X1  (x). 

In  order  to  account  for  the  presence  of  the  crack,  the  available  mode  shapes  of  the 
cracked  beam  are  not  applicable  as  they  are  not  twice  differentiable.  Therefore,  a  new 
formulation  of  the  mode  shapes  of  the  cracked  beam  is  proposed  and  modeled  in  the  present 
work.  In  the  proposed  modeling  approach,  Rayleigh-Ritz  approximation  is  applied  in 
formulating  approximate,  twice  differentiable  mode  shapes  of  the  cracked  beam. 

In  the  Rayleigh-Ritz  method,  mode  shapes  of  the  vibrating  beam  are  approximated  as 

n 

j(fc)(x)  =  ^  c\k^(pl{pc)  (9) 

i= 1 

where cf^’s  are  unknown,  constant  Ritz  coefficients  and  <pl(x)  are  the  admissible  functions, 
satisfying  boundary  conditions  of  the  system.  In  the  proposed  approximation,  the  crack  is  still 
assumed  to  be  treated  as  a  massless  rotational  spring  with  the  resultant  compliance  of  Cc  = 

11,85/ife  ( ac/hb )2 .  When  using  the  mode  shapes  of  the  un-cracked  beam  as  the  admissible 

Eb^b 

functions  for  Rayleigh-Ritz  approximation,  the  maximum  strain  energy  needs  to  be  modified  to 
account  for  the  loss  of  energy  due  to  the  presence  of  the  crack.  This  energy  loss  is  proportional  to 
the  amount  of  the  added  flexibility  due  to  the  presence  of  crack/rotational-spring. 

For  the  beam  with  a  uniform  cross  sectionthe  modified  stiffness  and  mass  matrices  for  the 
transverse  vibrations  of  the  beam  after  considering  the  effect  of  presence  of  the  crack  are 
modified  as 


-i ij  =  [  EiJhtp1  (x)4>j  (x)dx  -  EbIbCE[(pl  (xc)][(pi  (xc)] 
Jo 


rrii 


-s' 

J0 


(10) 


p  hb  b  <p 1  (x)  (x)  dx 


Therefore  the  approximation  problem  will  be  reduced  to  finding  eigenvalues  and 
eigenvectors  of  the  following  matrix  equation 
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|[/c]  —  (to^^)2[m] 


=  o 


(11) 


wherec^l  =  {c[k\  c^\  ...  ,  j  is  the  vector  of  Ritz  coefficients.  Once  vector  is  found 

the  mode  shapes  of  the  cracked  beam,  X^k\x),  can  be  calculated. 

By  performing  numerical  analysis,  it  is  shown  that  the  proposed  mode  shape  formulation 
provides  a  close  approximation  of  the  mode  shapes  of  the  vibrating  cracked  beam  compared  to 
the  available  cracked  beam  modeling  approaches. 

For  each  depth  and  position  of  the  crack,  the  modeling  approach  proposed  in  the  present 
work  is  used  for  approximating  the  resonant  frequencies  of  the  free  vibrations  of  the  beam.  These 
frequencies  of  the  first  three  modes  of  vibrations  versus  the  crack  position  along  the  beam  are 
plotted  in  Figure  26  -  Figure  28  and  compared  to  the  available  cracked  beam  modeling 
approaches.  In  order  to  have  a  better  picture  of  the  results,  graphs  for  only  two  crack  depth  ratios 
of  ac/hb  =  0.2  and  ac/hb  =  0.4  are  shown  in  these  figures.  In  these  graphs:  a)  long  dash  line 
represents  frequencies  of  the  reference  (available)  modeling  approach  for  crack  depth  ratio  of 
ac/hb  =0.2,  b)  dash  dot  line  represents  frequencies  of  the  proposed  modeling  approach  for 
crack  depth  ratio  of  ac/hb  =  0.2,  c)  dash  line  represents  frequencies  of  the  reference  modeling 
approach  for  crack  depth  ratio  of  ac/hb  =  0.4  and  d)  solid  line  represents  frequencies  of  the 
proposed  modeling  approach  for  crack  depth  ratio  of  ac/hb  =  0.4. 


Figure  26.  Resonant  frequency  of  the  1st  mode  of  the  beam  vs.  crack  position  along  the  beam  for 
a)  simply  supported  and  b)  fixed-fixed  boundary  conditions. 
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Figure  27.  Resonant  frequency  of  the  2nd  mode  of  the  beam  vs.  crack  position  along  the  beam  for 
a)  simply  supported  and  b)  fixed-fixed  boundary  conditions. 


a) 


Crack  Position  (xc.  L) 


Figure  28.  Resonant  frequency  of  the  3rd  mode  of  the  beam  vs.  crack  position  along  the  beam  for 
a)  simply  supported  and  b)  fixed-fixed  boundary  conditions. 


Fatigue  damage  characterization  using  laser  ultrasonic  methods 

Laser  ultrasonics  deals  with  the  generation  and  detection  of  ultrasound  in  a  solid,  liquid, 
or  gaseous  medium  using  a  laser  light.  This  type  of  system  is  particularly  attractive  to 
nondestructive  structural  and  materials  characterization  of  solids  due  to  the  nature  of  non- 
contact,  non-destructive,  and  in-situ  measurements.  A  laser  ultrasonic  based  SHM  method  has 
been  developed  in  this  project.  To  isolate  shear  stiffness  anisotropy  variations,  entirely  new  laser 
ultrasonic  methods  based  on  laser  line  sources  to  provide  linearly  polarized  shear  waves  have 
been  developed.  These  methods  have  been  used  to  assess  material  anisotropy  using  shear  wave 
birefringence  measurements  in  rolled  materials  and  can  be  used  to  quantitatively  isolate  various 
orientation  distribution  coefficients.  The  simplicity  of  this  approach  and  its  robust  character 
provide  entirely  new  directions  for  ultrasonic  sensing  of  fatigue  damage  in  metal  alloys. 
Successful  implementation  of  this  experimental  technique  has  been  used  to  measure  inherent 
material  anisotropy  in  aluminum  and  copper  alloys. 
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Experimental  equipment,  set  up  and  procedures 

Some  of  the  waveform  variations  recorded  with  the  previous  testing  methods  used  in  the 
program  on  lug  joints  indicated  that  changes  were  occurring  in  the  material’s  elastic  properties  in 
the  fatigue  damaged  region.  Based  on  analysis  of  a  subset  of  these  data  sets,  we  decided  to 
investigate  shear  anisotropy  variations  that  might  be  related  to  fatigue  damage  processes.  To 
isolate  these  variations,  we  developed  new  methods  using  linearly  polarized  shear  waves  to 
measure  wave  speed  variations  as  a  function  of  material  anisotropy.  However,  since  most  rolled 
plate  materials  exhibit  a  baseline  anisotropy  related  to  preferred  grain  orientation,  a  preliminary 
study  was  performed  to  characterize  shear  birefringence  effects  related  to  rolling-induced 
anisotropy. 
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Figure  29.  Ultrasound  generated  by  laser  line  source  with  contacting  transducer  for  detection. 
Orienation  of  sample  relative  to  the  line  source  is  also  shown. 


A  pulsed  laser  was  used  for  ultrasonic  generation  with  the  beam  focused  into  a  narrow 
line.  To  create  the  line,  a  double  concave  lens  and  double  convex  lens  were  used  to  expand  and 
collimate  the  beam,  while  a  cylindrical  lens  focused  the  beam  to  a  line  (Figure  29).  The  line 
source  forces  a  linear  polarization  of  the  shear  wave  which  is  perpendicular  to  the  orientation  of 
the  line.  The  line  source  is  centered  on  the  sample  which  is  mounted  to  a  rotational  stage  such 
that  0  and  180  degrees  are  interchangeable.  Using  the  experimental  geometry  shown  in  Figure 
29,  all  ultrasonic  modes  have  a  common  direction  of  propagation  such  that  rotation  of  the  sample 
does  not  affect  the  longitudinal  wave  since  its  polarization  and  propagation  directions  are 
unchanged  (both  are  in  the  through-thickness  direction  of  the  specimen).  In  the  first  series  of 
measurements  using  this  experimental  approach,  the  conical  transducer  was  used  for  detection  to 
quickly  characterize  the  system  and  generate  data.  The  transducer  was  mounted  on  the  back  side 
of  the  specimen  using  a  gel  couplant  so  that  the  aperture  area  was  centered  and  the  sample  could 
be  easily  rotated  over  the  transducer  while  still  making  sufficient  contact.  Data  were  collected  at 
15  degree  intervals  on  brass  and  copper  (2”x2”x0.1875”  square,  rolled  plates)  as  well  as  an 
aluminum  alloy  (2”x2”x0.25”  square,  rolled  plates).  Although  brass  and  copper  were  not  the 
focus  of  this  study,  they  were  added  to  the  sample  set  as  materials  with  known  texture  and 
detectable  anisotropic  characteristics.  Initial  measurements  indicated  that  rolling-related 
anisotropy  could  be  sensed  but  quantitative  measurements  were  not  obtained.  Consequently,  it 
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was  necessary  to  implement  interferometric  detection  (Figure  30)  to  provide  highly  repeatable, 
high  fidelity  measurements. 
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Figure  30.  Laser  line  source  for  ultrasonic  generation  and  interferometer  for  measuring 
displacement,  shown  with  the  orientation  of  the  sample  relative  to  the  line  source. 


Laser  line  source  experiments 
using  the  interferometer  required  that  the 
specimen  be  mounted  on  the  rotational 
stage  so  that  the  specimen  faces  were 
perfectly  perpendicular  to  the 
interferometric  beam  to  minimize  the 
angle  of  reflectance  from  the  incident 
beam  path.  The  interferometric  beam 
was  focused  on  the  opposite  side  of  the 
specimen  from  the  laser  source  and  was 
centered  so  that  interrogation  of  the 
region  remained  constant  as  the 
specimen  rotated.  The  same  data 
acquisition  methods  were  used  that  were 
described  for  the  transducer. 

Laser  line  source  results 


Plots  of  shear  wave  speed  as  a  function 
of  sample  orientation  for  each  sample  are 
shown  in  Figure  31  and  Figure  32. 

Brass  and  copper  samples  were 
tested  to  provide  insight  to  the  operation 
of  the  laser  line  source  experiment. 


Figure  31.  Series  of  copper  alloys.  A  double  y-axis 
plot  is  used  to  show  extruded  brass  velocities  are 
faster  than  the  rolled  plates  of  copper  and  brass. 


Angle  of  polarization  to  rolling  direction 

Figure  32.  Series  of  aluminum  alloys  from  laser  line 
source-in/laser-out  method.  A12024  exhibits 
measureable  anisotropy,  while  6061  appears  to 
maintain  isotropic  behavior. 


ave  speeds  c 
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Brass  is  known  to  exhibit  high,  rolling-induced  anisotropy  which  is  well-documented  in  the 
literature  and  can  be  used  to  establish  the  ability  of  this  technique  to  assess  anisotropy.  For 
brass,  a  greater  than  4.5%  difference  in  shear  wavespeed  between  maximum  and  minimum  was 
measured  with  an  uncertainty  of  +/-0.3%. 

Figure  32  displays  velocity  results  for  aluminum  alloys  2024  and  6061.  For  A12024  (in 
both  rolled  plate  and  extruded  bar  form)  maximum  velocities  can  be  clearly  noted  at  0  and  1 80 
degrees  from  the  rolling  direct,  while  A16061  is  isotropic  with  respect  to  shear  wave  polarization. 
Overall,  uncertainty  for  extracting  shear  wave  speeds  was  +/-0.3%  while  the  maximum  variation 
in  wave  speed  for  A12024  was  1.5%  -  significantly  above  the  measurement  uncertainty. 

Task  3:  Fatigue  Life  Prediction  Using  Hybrid  Prognosis 

Current  practice  in  fatigue  life  prediction  is  based  on  assumed  initial  structural  flaws 
regardless  of  whether  these  assumed  flaws  actually  occur  in  service.  Small  deviations  of  the 
initial  conditions  and  model  parameters  may  generate  large  errors  in  the  expected  dynamical 
behavior  of  fatigue  damage  growth.  This  project  develops  an  integrated  approach  of  SHM  and 
adaptive  prognosis  model  that  not  only  estimates  the  current  health,  but  can  also  forecast  the 
future  health  and  calculate  RULE  of  an  aerospace  structural  component  with  high  level  of 
confidence. 

Generalized  Bayesian  data  driven  model  for  damage  prediction 

The  lack  of  general 
applicability  of  available  physics 
based  crack  growth  model  is  due  to 
the  complex  dependency  of  crack 
propagation  on  different  factors  such 
as  variability  due  to  (i)  material  (e.g., 
microstructure)  (ii)  geometric 
complexity  (iii)  manufacturing 
process  (e.g.,  heat  treatment,  cold 
deformation)  (iv)  loading  (v) 
environmental  effects  (e.g., 
temperature  and  humidity).  To 
incorporate  these  factors  in  a  damage 
propagation  model,  a  generalized  Bayesian  probabilistic  framework  can  be  employed.  The 
Bayesian  framework  is  a  data  driven  approach.  The  schematic  of  a  probabilistic  Bayesian 
framework  for  future  step  damage  state  forecasting  is  shown  in  Figure  33.  The  goal  of  a 
probabilistic  Bayesian  forecasting  approach  is  to  compute  the  posterior  distribution  of  an 
unknown  target,  i.e.,  to  predict  f(yn+ ±\D  =  {Xi,y;}i=o where  yn+1  represents  either 
the  crack  length  an+1  or  the  crack  growth  rate  (da/dN)n+l  at  damage  level  n  +  1.  The  (n  +  l)th 
damage  level  distribution  is  obtained  based  on  available  input-output  data  up  to  the  current 
damage  level  and  future  input  parameter  xn+l  .  It  is  noted  that  the  input  Xj  at  any  fatigue  damage 
level  i  is  a  vector  and  constitutes  the  various  fatigue  affecting  parameters.  In  the  Bayesian 
framework,  the  predictive  distribution  of  a  target  can  be  found  by  conditioning  the  targets 

yo, . ,yn-i,yn,  yn+i  that  are  affected  by  the  corresponding  random  inputs  xo, . ,xn-i,xn,  xn+i. 

Now  we  can  define  a  priori  over  the  space  of  possible  functions  to  model  the  target  (crack  length 


Figure  33.  General  Bayesian  probabilistic  framework 
for  future  step  damage  state  forecasting 
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or  crack  growth  rate)  as/(y|a),  where  a  represents  parameters  that  account  of  modeling 
uncertainties  in  the  form  of  curve  fitting.  It  is  assumed  that  modeling  the  uncertainty  parameters 
a  can  account  for  the  effect  of  loading  interaction  (e.g.,  retardation  effect)  in  addition  to 
modeling  uncertainties.  A  priori  noise  function  /($|/?)  can  also  be  defined,  where  9  is  some 
appropriate  noise  function  that  arises  due  to  scatter  in  material  micro  structure  and  P  is  another 
set  of  hyperparameters  used  to  model  the  uncertainty  due  to  scatter.  Now  if  the  parameters  a  and 
P  are  known,  the  conditional  probability  can  be  expressed  as 

f{yn\[Xi=i . ,n,tX'P})  =  J(yn|*s=i. . ,n>y'V)f(y\a)fW)dydd  (12) 

where  yn  =  {yo,  . ,  yn-i,  yn}  and  0  denotes  the  underlying  function  which  corresponds  to  the 

target  functions  and  noise  due  to  microstructural  scattering  respectively.  Since  yo, . ,yn-i,  and  yn 

are  conditioned  random  variables  in  the  observed  set  of  targets,  the  conditional  distribution  of 
yn+i  can  be  written  as  follows 

^  f  ( -» 7  .  ^  \  -V  .  ^  ^  A  7  -l9  ^ 

(13) 


/(y„«|0  =  fey<}i=i . 


Damage  prediction  using  data  driven  Gaussian  Process  (GP)  approach 

The  expression  given  in  Eq.  13  is  accurate 
if  there  exists  a  ’’linear  dependence”  between  the 
two  random  variables  Xj  and  xj.  However,  in 
general  for  fatigue  damage  modeling,  the  input 
variables  xi  and  xj  may  not  have  a  linear 
relationship.  The  nonlinearity  is  more  pronounced 
during  the  unstable  and  transient  damage  growth 
regime.  As  parameters,  such  as  crack  length  and 
loading  information,  are  introduced  in  the  input 
space,  the  relation  between  input  variables  x;  and 
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Xj  becomes  more  complex.  To  avoid  the 


Figure  34. 

dimensional 

space. 


Schematic  of  high- 
transformation  of  input 


nonlinearity  problem,  Gaussian  process  (GP)  is 

used  to  map  the  nonlinear  parameters  into  a  high  dimensional  space,  as  shown  in  Figure  34.  It  is 
parameterized  by  a  mean  and  a  covariance  function.  With  state  information  available  up  to  nth 


damage  level,  the  predictive  distribution  at  (n+l)th  damage  level  can  be  given  as: 


P(an+i,an)  =  -  expljp- 


(an  + 1  Un+l)1 


2  a 


■) 


an  + 1 


an+1  =  kTK-x' 


an  + 1 


=  k  -  kTK-J 


(14) 

(15) 

(16) 


where  an+1  is  the  predictive  mean  at  the  (n+l)th  damage  level,  and  agn+1  is  the  associated  error 

variance  of  the  prediction,  Kn  is  the  n  x  n  kernel  matrix  for  the  vector  an ,  k  and  kT  are  the 
partitioned  components  of  kernel  matrix. 

For  the  future  damage  state  prediction,  the  GP  model  given  by  Eq.  14-16  recursively 
predicts  the  future  damage  states  based  on  the  last  on-line  data  available.  The  training  data  set  D 
and  test  input  vector  an+„  can  be  written  as, 
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D  = 


- 

training 

data  matrix 

target  vector  - 

ao 

ai 

a2 

ad— 1 

ad 

ai 

a2 

a3 

"•  ad 

ad+l 

an— d 

an— d+1 

a3n— d+2 

an— 1 

an 

an  — d+1 

an-d+2 

a3n— d+3 

•"  an 

qP 
n  +  l 

an— d+1 

an— d+2 

a3n— d+3 

. . .  a? 

n  +  l 

ap 

dn+2 

P  P 

_dn— d— 1+n  dn— d+n 

ap 

n-d+l+n 

...  aP 

n— 2+fT 

ap 

dn-l+n  -1 

test  input  data  vector 

an+n  — 

ap 

dn— d— l+n 

ap  ap 

an-d+n  3n- 

-d+l+n 

ap 

n— 2+n“ 

(17) 


(18) 


where  in  Eqs.17  and  18,  the  subscript  n  is  the  damage  level  to  which  the  last  on-line  data  is 
available,  the  subscript  n  is  the  damage  level  after  the  availability  of  the  last  on-line  data  and  the 
superscript  p  indicates  the  predicted  damage  index. 

Hybrid  prognosis  model 

Fatigue  crack  growth  behavior  for  a  given  specimen  can  be  predicted  by  combining 
knowledge  of  the  underlying  mechanics  of  crack  growth  and  future  loading.  The  hybrid 
prognosis  framework  presented  in  this  paper  considers  simple  crack  growth  models  whose 
behaviors  are  inferred  and  updated  using  data-driven  approaches.  The  combination  of  physics 
and  data-driven  approaches  allows  for  the  consideration  of  proper  fracture  mechanisms  while 
correcting  for  material  variations  and  uncertainty  in  the  model  parameters  using  data-driven 
model  updating.  Thus,  although  simple  physics  models  are  used,  the  accuracy  of  the  hybrid 
framework  is  greater  than  those  of  data-driven  and  physics  based  models  alone,  as  shown  in  the 
results  presented  in  a  later  section. 

Linear  elastic  fracture  mechanics  and  most  fracture  theories  state  that  the  crack  growth 
rate  ( dal dN )  is  a  function  of  the  stress  intensity  factor  range  ( A K ),  as  shown  in  Eq.  19. 

—  =  f(AK)  (19) 

dF  V  ’ 

where,  K  =  K  max  K  min 

Due  to  the  exponential  nature  of 
crack  growth,  most  models  typically 
describe  the  relationship  between  crack 
growth  rate  and  stress  intensity  factor 
(SIF)  using  log-log  transforms.  The 
commonly  observed  trend  showing 
three  critical  zones,  stages  I-III  of  crack 
growth  is  shown  in  Figure  35. 

Prediction  initiation  in  stage  I  can  often 
result  in  large  errors  with  respect  to  life 
since  cracks  can  grow  on  the  order  of 
N3  or  N4.  Typically,  prognosis 
algorithms  are  applied  during  stage  II 


Stress  intesity  factor  range,  AX,  (MPaVm),  log  scale 

Figure  35.  The  relationship  between  crack  growth 
rate  and  SIF  is  highly  nonlinear  even  after  log-log 
transform.  There  is  a  stable  linear  portion,  indicated 
by  the  black  line  that  is  representative  of  Stage  II  or 
sub-critical  crack  growth. 
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or  sub-critical  crack  growth  and  are  used  to  predict  ultimate  fracture.  In  some  cases,  for 
example,  constant  amplitude  loading,  this  regime  is  linear.  Models  such  as  Paris’  Law  are  well 
suited  to  capture  this  behavior.  For  cases  such  as  overloads  and  under-loads,  this  regime  can  be 
highly  nonlinear  and  discontinuous  requiring  the  use  of  advanced  physics  models,  which  are 
often  unavailable.  In  the  hybrid  prognosis  framework  presented  here,  the  exact  relationship 
between  crack  growth  rate  and  SIF  for  a  given  cycle  is  inferred  from  the  available  data  based  on 
an  assumption  of  a  linear  relationship  with  non-constant  coefficients  in  log  space.  The  constants 
of  the  linear  fit  are  a  function  of  historical  crack  growth  data,  future  loading  (i.e., 
overloads/under-loads),  and  basic  material  properties  and  cycles,  continuously  evolving  and 
adapting,  as  more  data  are  available.  Initial  estimates  for  these  parameters  (i.e.,  prior  to  data 
acquisition)  can  be  obtained  through  the  basic  material  constants  used  in  Paris’  Law,  which 
reduces  the  crack  growth  rate  estimation  to  a  classical  Paris’  Law  extrapolation.  However,  since 
the  SHM  framework  provides  data  on  crack  length  and  locations,  as  well  as  load  monitoring  and 
cycle  counting,  the  constants  are  updated,  allowing  them  to  model  and  capture  the  nonlinear  and 
discontinuous  behavior.  In  order  to  predict  the  fatigue  crack  growth  of  a  specimen,  Eq.  19  needs 
to  be  formulated  in  terms  of  measureable  parameters  and  integrated  until  ultimate  fracture. 
Therefore,  the  parameters  in  Eq.  19  must  be  written  in  terms  of  these  data,  requiring  that  SEF  is 
related  to  known  or  quantifiable  parameters.  We  assume  that  the  SIF  can  be  expressed  as  a 
general  function 

A KN=f(aN,PN,S)  (20) 

where,  S  is  a  geometric  parameter.  For  simple  structures,  analytical  expressions  of  SIF  are 
available  that  describe  its  dependence  on  geometry,  crack  length,  and  applied  load.  When  an 
analytical  expression  is  available,  it  can  then  be  directly  substituted  into  Eq.  20  and  the  future 
crack  growth  can  be  calculated.  However,  in  the  absence  of  this  information,  numerical  methods 
must  be  utilized  to  provide  estimates  of  SIF.  Either  method  is  acceptable  and  suitable  for  use  in 
the  proposed  framework.  The  only  model  parameter  necessary  to  apply  the  proposed  prognosis 
model  is  the  stress  intensity  factor  (SIF).  A  SIF  mapping  can  be  created  for  any  specimen  as  a 
function  of  the  crack  tip  location.  Finite  element  simulations  can  be  run  for  different  crack  tip 
locations,  which  will,  in  turn,  serve  as  the  training  data.  A  regression  model  can  then  be  used  to 
map  the  data,  allowing  for  evaluation  of  the  SIF  for  any  given  crack  tip  location  obtained  from 
the  experiments. 


Experimental  validation 

Validation  using  cruciform  testing  articles 

To  numerically  validate  the 
integrated  prognosis  algorithm,  a  fatigue  test 
was  performed  on  an  Al-6061  cruciform 
specimen  under  biaxial  loading.  The  loaded 
cruciform  specimen  in  an  MTS  biaxial 
fatigue  test  frame  can  be  seen  in  Figure  36. 
The  specimen  was  subjected  to  constant 
amplitude  fatigue  loading  with  maximum 
amplitude  (s max)  4  kips  and  load  ratio  of  R 
=  0.1,  and  the  biaxial  machine  actuator  was 
operated  with  a  frequency  of  10  Hz.  It 
should  be  noted  that,  the  maximum  stress 


Figure  36.  Al-6061  cruciform  specimen  loaded 
in  a  MTS  biaxial  fatigue  test  frame. 


Cruciform 
specimen  under 
biaxial  loading 
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amplitude  was  equal  to  two  thirds  the  yield  stress  sF .  Based  on  nonlinear  finite  element  analysis 
of  cruciform  specimen,  the  yield  stress  was  approximated  as  sF  =  6kips.  Also  note  that  both  the 
x-axis  actuator  and  y-axis  actuator  of  the  biaxial  frame  were  subjected  to  in-phase  fatigue 
loading.  For  on-line  state  estimation,  passive  strain  gauge  sensors  were  used.  Two  strain  gauges 
were  mounted  on  the  web  area,  one  strain  gauges  mounted  on  the  horizontal  flange,  and  the  other 
one  on  the  vertical  flange  of  the  cruciform  specimen.  In  addition,  a  hole  in  the  center  of  the 
specimen  was  made  to  create  crack  initiation  in  the  web  area  of  cruciform  specimen.  To 
accelerate  damage  growth  an  EDM  notch  of  1  mm  length  was  made  at  left  bottom  quadrant 
boundary  of  the  central  hole  (45 o  to  the  vertical  axis).  A  48  channel  NI  PXI  system  was  used  to 
collect  the  strain  gauge  signals  and  the  measurements  from  the  biaxial  machine  load  cells.  In 
addition,  a  high  resolution  SONY  camera  was  used  to  visually  monitor  the  crack  growth.  The 
data  acquisition  system  and  the  computer  capturing  the  visual  image  were  synchronized  with  the 
biaxial  machine  controller  to  collect  the  time  synchronized  data /  image  at  a  specified  interval  of 
DA  =  1500 cycles.  The  data  and  image  collection  started  at  approximately  Ilk  cycles.  The  image 
and  sensor  data  were  collected  at  47 
different  time  instances.  For  the  first 
44  instances,  the  signals  and  images 
were  collected  while  the  biaxial 
machine  was  running  and  during  the 
last  three  instances  the  data  was 
collected  when  the  machine  was  not 
running.  This  leads  to  a  total  of  44 
different  damage  cases  with  the  last 
damage  state  occurring  at  75.5  kcycles. 

Figure  37  shows  the 
comparison  between  single  step  ahead 
forecasted  state  and  actual  damage 
state  (or  damage  index)  with  on-line 
data  available  up 
to  the  previous 
damage  level.  As 
seen  in  the  figure, 
the  prognosis 

algorithm  starts 
predicting  from 
estimated  from 
the  sensor  signals. 

The  threshold 

value  of  0.7  is  70% 
of  the  final 
damage  index 
value  of  one.  The 
damage  index 
reaches  its  final 

value  of  one  Figure  38.  Multi-step  ahead  damage  state  prediction  using  off-line 
when  there  is  no  predictive  model. 


Figure  37.  One-step  ahead  damage  state  prediction 
using  off-line  predictive  model. 
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cross-correlation  between  the  input  u  and  y.  This  is  because  the  specimen  undergone  complete 
failure.  It  should  be  noted  that  choosing  the  critical  damage  index  value  of  0.7  was  based  on  the 
results  from  previously  performed  similar  experiments. 

Unlike  the  single  step  ahead  prediction,  the  multi  step  ahead  prediction  recursively 
predicts  the  damage  state  multiple  steps  ahead  of  the  damage  level  at  which  last  online  data  was 
available.  Figure  38  shows  the  multi  step  ahead  state  prediction.  Similar  to  single  step  ahead 
prediction  process,  the  prognosis  algorithm  was  started  after  the  6th  damage  level  (i.e.  at  18.5 
kcycles).  From  the  7th  damage  level  (from  20  kcycles),  damage  indices  were  predicted  and  then 
fed  back  to  the  prognosis  model  to  update  the  Gaussian  process  training  data  matrix  and  the  test 
input  vector.  The  feedback  process  and  the  corresponding  future  state  predictions  were  continued 
recursively  as  long  as  the  predicted  damage  index  did  not  reach  its  critical  value  of  0.7.  It  is  to  be 
noted  that  unlike  the  single  step  ahead  prediction,  the  training  data  matrix  and  the  corresponding 
test  input  vector  were  updated  with  off-line  model  predicted  states,  rather  than  being  updated 
with  on-line  model  estimated  states,  which  could  not  be  available  in  real  time.  It  can  be  seen 
from  Figure  38  that,  with  on-line  data  available  up  to  damage  level  23  (at  44  kcycles),  the  multi 
step  ahead  predicted  states  fails  to  reach  the  critical  value  of  0.7.  This  is  because  the  predictive 
model  was  unable  to  learn  the  damage  growth  dynamics.  It  is  also  to  be  noted  that  if  the 
predictive  model  does  not  learn  the  damage  growth  dynamics  it  keeps  on  running  with  only 
predicting  unvarying  damage  indices.  The  predicted  unvarying  damage  indices  time  series  can 
also  be  seen  from  the  Figure  38.  Without  satisfying  the  threshold  criteria,  the  prediction  of 
unvarying  damage  indices  could  have  continued  indefinitely.  However  to  reduce  the 
computational  expenses,  the  prognosis  algorithm  was  stopped  at  certain  times.  The  criteria  for 
stopping  the  algorithm  was  if  the  rate  of  damage  index  growth  was  not  greater  than  1x10' 
1=cycles  for  six  consecutive  damage  levels,  the  off-line  predictive  model  had  to  be  terminated. 
This  was  because  of  physical  reason,  if  the  damage  growth  was  slow  enough,  the  predicted 
damage  index,  would  never  reach  the  critical  value  even  if  the  algorithm  had  to  run  indefinitely. 
From  Figure  10  it  is  also  seen  that,  the  first  multi  step  ahead  prediction  curve,  that  reaches  the 
critical  value  starts  from  damage  level  24  (from  45.5  kcycles).  Beyond  this  damage  level,  the 
multiple  step  ahead  prediction  increasingly  converges  with  the  actual  damage  index.  From  the 
above  mentioned  observations,  it  can  be  assumed  that  the  prediction  horizon  (or  the  true  positive 
regime)  was  between  damage  level  24  (45.5  kcycles)  and  damage  level  42  (72.5  kcycles),  during 
which,  the  predicted  damage 
states  reached  its  critical 
value. 

The  residual  useful 
life  at  any  given  damage 
level  (up  to  which  the  last 
online  data  was  available) 
was  estimated.  Figure  1 1 
shows  the  comparison  of 
predicted  RULE  and  actual 
RULE.  From  the  figure  it 
can  be  seen  that,  there  is  a 
good  correlation  between 
predicted  and  actual  RULE 

in  the  true  positive  regime  Figure  39.  Comparison  of  predicted  RULE  and  actual  RULE. 
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i.e.,  between  45.5  kcycles  and  72.5  kcycles.  Also  as  more  and  more  online  data  becomes 
available,  better  correlation  between  predicted  RULE  and  actual  RULE  is  observed. 

Validation  using  lug  joint  test  articles 

For  the  validation  of  the  method  on  a  real-time  specimen,  an  Aluminum  2024-T351  lug 
joint  subjected  to  fatigue  loading  was  instrumented  and  interrogated.  The  experimental  setup  is 
the  same  as  reported  in  last  two  sections.  To  track  the  crack  growth,  two  cameras  were  mounted 
on  to  the  frame,  each  focusing  on  the  crack  on  the  front  and  rear  of  the  specimen.  The  captured 
images  were  used  to  calculate  the  crack  length. 

In  order  to  evaluate  the  stress  intensity  factor  (SIF),  a  finite  element  simulation  was  run 
using  Abaqus.  A  full  3D  finite  element  model  of  the  lug  joint  was  created  in  Abaqus.  The  crack 
propagation  direction  was  considered  as  normal  to  the  crack  front  plane  for  the  SIF  calculation. 


Figure  40.  Finite  element  model  of  the  lug  joint  with  crack 

A  grid  of  15mm  x  25mm  was  made  at  one  shoulder  of  the  lug  joint,  and  the  crack  tip  was 
modeled  for  17  different  locations  on  the  grid.  The  left  pinhole  on  the  lug  joint  was  fixed  in  all 
directions,  and  the  right  pinhole  was  allowed  to  move  along  the  direction  of  the  loading.  The 
SIF  was  calculated  with  a  load  of  13kN.  It  was  also  observed  that  the  SIF  varies  linearly  with 
loading.  It  is  noted  that  for  the  first  crack  tip  location  (109.00,  79.50),  as  the  load  reduces  10 
times,  the  SIF  also  reduces  to  l/10th  of  the  original  value.  For  the  crack  tip  location  (104.00, 
79.50),  as  the  load  reduces  by  Vi,  the  SIF  also  reduces  by  Vi.  And,  for  the  crack  tip  location 
(113.16,  84.53),  as  the  load  reduces  by  l/4th,  the  SIF  also  reduces  by  l/4th.  These  observations 
indicate  that  the  SIF  varies  linearly  with  the  load.  Hence,  for  any  given  loading,  the  SIF  can  be 
calculated  as  a  linear  function  of  the  load.  In  our  experiments,  the  SIF  (K)  was  calculated  for  unit 
load.  Then,  for  any  given  loading,  the  calculated  SIF  was  multiplied  by  the  load  to  get  the  new 
SIF  for  that  particular  load.  From  this  process,  we  are  able  to  find  the  SIF  as  a  function  of  the 
crack  tip  location  for  different  loads.  After  the  SIF  was  evaluated,  a  mapping  for  the  SIF  was 
created  as  a  function  of  the  crack  tip  location.  Once  the  mapping  was  created,  the  SIF  can  be 
evaluated  from  the  mapping  for  any  obtained  crack  tip  location  during  the  experiment.  To  create 
a  mapping  of  the  SIF,  three  methods  were  used.  First,  the  mapping  was  created  using  the  surface 
fitting  toolbox  in  MATLAB,  and  then  using  LASSO,  and  then  using  RVM.  The  first  two 
methods  are  deterministic,  whereas  RVM  is  a  probabilistic  method,  which  provides  the 
confidence  intervals  of  the  prediction. 

SIF  mapping  using  RVM 

For  the  SIF  mapping  using  RVM,  the  input  parameters  are  the  crack  tip  coordinates  (x, 
y),  and  the  output  parameter  is  the  SIF  for  a  given  crack  tip  location.  While  the  LASSO  is  a 
useful  tool  for  robust  regression,  it  yields  a  point  estimate  for  the  regressed  model.  RVM  is  a 
powerful  technique  that  can  perform  full  probabilistic  regression,  and  it  has  the  advantage  of 
additionally  providing  a  measure  of  uncertainty  (confidence  intervals)  on  the  regressed  estimate. 
Thus,  using  the  RVM  for  the  SIF  mapping,  a  probability  distribution  on  the  SIF  as  a  function  of 
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x  and  y  position  can  be 
obtained.  In  this  setting, 
the  uncertainty  in  SIF  can 
be  transferred  naturally  to 
uncertainty  in  crack 
growth  rate,  for 
subsequent  use  in  the 
prognosis  model.  In  our 
evaluation,  the  inputs 
used  included  17  data 
points,  and  the 
predictions  were  plotted 
on  a  grid  of 
201x301=60501  data  points.  For  each  of  the  60501  grid  points,  a  lower  bound  and  an  upper 
bound  were  established,  which  can  be  related  to  the  variance  of  the  prediction.  The  input  data 
was  mapped  with  RVM  using  a  Gaussian  kernel,  and  the  obtained  fit,  with  lower  bound  and 
upper  bound  surfaces,  is  shown  in  Figure  41.  The  surfaces  were  plotted  with  a  variance  of  2a. 
For  this  set  of  input  and  output  data,  Gaussian  kernel  was  observed  as  providing  good  results; 
hence  this  kernel  was  used  for  the  mapping. 


_ 25 

E 
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Figure  41.  Plot  of  surface  fit  using  RVM. 
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